/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

///\file
/// Definitions of particle classes. \ingroup makegrid

// $Id$

#ifndef __sphparticle__
#define __sphparticle__

#include "constants.h"
#include "mcrx-types.h"
#include "blitz/array.h"
#include <string> 

namespace mcrx {
  template <typename cell_data_type > class grid_cell;
  class particle_base;
  template <typename data_type> class particle;
}

template< typename, typename, int> class interpolator;

// *** class particle_base ***

/** Base class for particles.  This class knows about how to calculate
    the overlap between particles and grid cells.  It's a base class
    for the concrete particle classes, depending on their
    content.  \ingroup makegrid */
class mcrx::particle_base {
public:
  typedef mcrx::T_float T_float;
  typedef mcrx::vec3d vec3d;
private:
  int idnum; ///< Particle ID number.
  vec3d pos; ///< Particle position.
  vec3d vel; ///< Particle velocity.
  T_float h; ///< Particle radius (SPH smoothing length).

  //static ndim_interpolator* interpol;
  typedef interpolator < T_float, T_float, 4> T_interpolator;
  /** The interpolator used for overlap calculations.  Because the
      calculation of how much of a particle is inside a cube is an
      expensive three-dimensional integration, this is pre-calculated
      and instead an interpolation in a table is used.  */
  static T_interpolator* interpol;
  void create_interpolator ();
  
  T_float kernel (T_float) const;
  
public:
  particle_base (int i, const vec3d& p, const vec3d& v, T_float sz); 

  int id () const {return idnum;};
  const vec3d& position () const {return pos;};
  const vec3d& velocity () const {return vel;};
  T_float radius () const {return h;};
  T_float project (const vec3d& cellmin,const vec3d& cellmax) const;
  /// Returns true if the particle overlaps with the cubic region specified.
  bool overlap (const vec3d& cellmin,const vec3d& cellmax) const;
  /// Returns true if the particle center is within the cubic region specified.
  bool inside (const vec3d& cellmin, const vec3d& cellmax) const {
    return
      (cellmin [0] < position() [0]) && (cellmin [1] < position() [1]) &&
      (cellmin [2] < position() [2]) && (cellmax [0] > position() [0]) &&
      (cellmax [1] > position() [1]) && (cellmax [2] > position() [2]);};
  T_float draw_radius (T_float r) const;
  static std::string inp_dir;

  /** Calculates the kernel-weighted fraction of the particle inside
      the grid cell.  */
  template <typename cell_data_type>
  T_float project (const grid_cell<cell_data_type>& c) const {
    assert (all((abs(c.getmax()- c.getmin())/c.getsize()-1) < 1e-5));
    return project (c.getmin(), c.getmax());}

  /// Returns true if any part of the particle is inside the grid cell.
  template <typename cell_data_type>
  bool overlap (const grid_cell<cell_data_type>& c) const {
    return overlap (c.getmin(), c.getmax());}
  template <typename cell_data_type>
  bool inside (const grid_cell<cell_data_type>& c) const {
    return inside (c.getmin(), c.getmax());}

};

// *** class particle ***

/** A particle is a particle_base with some data attached.
    \ingroup makegrid */
template <typename data_type>
class mcrx::particle: public particle_base {
public:
  typedef data_type T_data;
private:
  T_data d; ///< The "content"  of the particle.

public:
  particle (int i, const vec3d& p, const vec3d& v, T_float sz, const T_data& dd);
  /// Returns the particle content.
  const T_data& data () const {return d;};
};


// ** constructors ***

/** Creates a particle with specified ID number, position, size and
    content.  */
template <typename data_type>
mcrx::particle<data_type>::particle (int i, const vec3d& p, const vec3d& v,
				     T_float sz, const T_data& dd ):
  particle_base (i, p, v, sz) , d (dd)
{}


/// Calculates the SPH smoothing kernel from Hernquist & Katz 1989.
inline
mcrx::particle_base::T_float
mcrx::particle_base::kernel(T_float r) const
{
  assert (r >= 0);

  const T_float c=1./(constants::pi*h*h*h);
  const T_float a=r/h;
  if(a<=1)
    return (1-1.5*a*a+0.75*a*a*a)*c;
  else if(a<=2)
    return (0.25*(2-a)*(2-a)*(2-a))*c;
  else
    return 0;
}



#endif
